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. Abstract 

Kolmogorov wave turbulence plays an important role for the 
CN ■ thermalization process following plasma instabilities in nonabelian 

gauge theories. We show that classical-statistical simulations in SU(2) 
Q>^ I gauge theory indicate a Kolmogorov scaling exponent known from 

CN ■ scalar models. In the range of validity of resummed perturbation 

theory this result is shown to agree with analytical estimates. We 
study the effect of classical-statistical versus quantum corrections and 
demonstrate that the latter lead to the absence of turbulence in the 



00 

■ far ultraviolet. 
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1 Introduction 



Heavy-ion collision experiments explore strongly interacting matter starting 
from a transient far-from-equilibrium initial state. Available data from the 
Relativistic Heavy Ion Collider reveals remarkable, unexpected properties 
such as rapid apparent thermalization with robust collective phenomena. 
The explanation of these findings from quantum chromo dynamics (QCD) 
provides a challenge for theory. There are theoretical indications that 
nonequilibrium dynamics related to plasma instabilities could be crucial for 
the understanding of the plasma's temporal evolution [H [21 [3], H]. These 
instabilities lead to exponential growth of occupation numbers in long 
wavelength modes on time scales much shorter than the asymptotic thermal 
equilibration time. After the fast initial period of exponential growth the 
dynamics slows down considerably. The system is still far from equilibrium 
at this stage and the subsequent slow evolution may be characterized by 
Kolmogorov wave turbulence 

In the turbulent range characteristic properties such as values for scaling 
exponents can be insensitive to the details of the underlying microscopic 
theory. The associated universality may lead to quantitative agreements for 
very different theories or energy scales. Wave turbulence has been studied 
in detail in weakly-coupled scalar theories in the context of early- universe 
reheating P, [7j. Though the underlying mechanisms are very different 
for QCD and for scalar inflaton dynamics, the subsequent evolution after 
an instability follows similar patterns: There a tachyonic or parametric 
resonance instability leads to the exponential growth of occupation numbers, 
followed again by a slow cascade P El 1101 [HI [12] • Wave turbulence in 
scalar theories is described by nonthermal distributions such as n{p) ~ 
p~'^ with K = 3/2 for cubic self-interactions in the presence of a field 
expectation value, or for quartic interactions k = 4/3 (5/3) describing 
particle (energy) cascades |6j. These distributions are distinct from a 
thermal high-temperature distribution ~ p~^. The nonthermal power- 
law behavior can be obtained from stationary solutions of perturbative 
Boltzmann equations neglecting quantum corrections [6j. A new class of 
nonperturbative scaling solutions, with strongly enhanced low-momentum 
fluctuations, has recently been found [7j. A very detailed picture of the 
different scaling solutions emerges using nonequilibrium renormalization 
group techniques [T3j. 

In contrast to the detailed picture for scalar quantum theories, our 
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knowledge about gauge theories is still in its infancies. Even at sufficiently 
high momenta, where perturbative estimates are expected to be applicable, 
conflicting statements can be found about the nature of power cascades 
following plasma instabilities [HI |T5] . It has been noted that the value for 
the exponents characterizing scaling behavior for gauge theories could deviate 
from those known for scalar fleld dynamics. For instance, a detailed numerical 
study in Ref. [H] for SU(2) gauge theory in Coulomb gauge suggests the value 
two for the scaling exponent, which is not fully understood perturbatively 
and seems to have no corresponding value in scalar models. The estimate 
employs a separation of scales between suitably deflned 'soft' and 'hard' 
momenta for sufficiently small characteristic running gauge coupling. The 
respective hard-loop effective theory of soft excitations is a nonabelian 
version of the linearized Vlasov equations of traditional plasma physics, 
which are based on coUisionless kinetic theory for hard particles coupled 
to a soft classical fleld p]. Vlasov equations can also be derived starting 
from classical-statistical fleld theory [16]. Therefore, they may also be 
considered as an approximation of the classical-statistical fleld theory limit of 
the respective quantum gauge theory [31 Sj. Classical-statistical lattice gauge 
theory provides a quantitative description in the presence of sufficiently large 
energy density or occupation numbers per mode. 

In this work we show that classical-statistical simulations following 
plasma instabilities in SU(2) gauge theory indicate a Kolmogorov scaling 
exponent known from scalar models (Sec. [2]). More precisely, we demonstrate 
that the numerical simulation is consistent with k = 4/3 characterizing 
particle cascades for 2^2 scattering associated to the phenomenon of 
weak wave turbulence. We show that this matches resummed perturbative 
estimates based on the two-particle irreducible (2PI) effective action (Sec. [3]). 
In this perturbative regime we study the effect of quantum corrections and 
demonstrate that turbulence is absent at sufficiently short distances. 

2 Classical-statistical lattice gauge theory 

Collisions of heavy nuclei leave behind a plasma of quarks and mostly gluons 
in a nearly flat region of space because of Lorentz contraction along the beam- 
or 2;-axis. If the quarks and gluons streamed freely into the surrounding space, 
then the distribution of particles would become locally highly anisotropic. In 
this case the stress tensor quickly acquires an oblate form with T^^ ~ Tyy ^ 
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T^^. The fastest growing mode of plasma instabilities then has its wave vector 
along the normal direction and generates a prolate contribution to the stress, 
which pushes the system towards greater isotropy IH |2l |3l H] . 

We describe this physics using classical-statistical Yang-Mills theory in 
Minkowski space-time following closely Ref . [1] , to which we refer for further 
technical details. It employs the Wilsonian lattice action for SU(N) gauge 
theory in Minkowski space-time: 

S[U] = -/?oEE{^(Trf^^,o. + Trf/i,o.)-l} 

X i<j ^ 

with X = x) and spatial Lorentz indices i,j = 1, 2, 3. It is given in terms 
of the plaquette variable 11^^^,^ = U^^fJJ^+f,^yUl_^^^^Ul^^, where Ul ,^^ = 11^^^,^ . 
Here Ux,fj. is the parallel transporter associated with the link from the 
neighboring lattice point x + fi to the point x in the direction of the lattice 
axis fi = 0,1,2,3. The definitions (3o = and /3, = 2Trl/(^2^) 

contain the lattice parameter 7 = ag/at, where denotes the spatial and at 
the temporal lattice spacings, and we will consider go = Qs = 9- 

We specify to SU(2) as the gauge group. The dynamics is solved in 
temporal axial gauge. Varying the action ([T]) w.r.t. the spatial link variables 
Ux,j yields the classical equations of motion. Variation w.r.t. to a temporal 
link gives the Gauss constraint. The coupling g can be scaled out of the 
equations of motion and we will set g = 1 for the simulations. We define the 
gauge fields a^ 

= -— Tr Kf/i(x)) (2) 

(Xg 

where cTq are the Pauli matrices with color index a = 1,2, 3. Correlation 
functions are obtained by repeated numerical integration of the classical 
lattice equations of motion and Monte Carlo sampling of initial conditions. 
To be specific, we consider the extreme anisotropy case described by an 
effectively 6{pz)-like initial Gaussian distribution. The initial time derivatives 
Af^{t = 0,x) are set to zero, which implements the Gauss constraint at all 
times. Shown results are from computations on A^^ = 128^ lattices, where 
gauge fixing to Coulomb gauge is employed using a stochastic overrelaxation 

-'^Here we use a different normalization of the gauge field than in equation (8) of Ref. [4j. 
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Figure 1: Fourier coefficients of the squared modulus of the gauge field in 
Coulomb gauge as a function of momentum for different snapshots in time. 



algorithm [T7] . We comment on the observed stability of the results on 64^ 
and 256^ lattices below. 

Fig. [1] shows the equal-time correlation function 

Fi'^^'^'it-p) = I d'xe-'^-{A';{t,^)Al{t,0)) (3) 

as a function of momentum for various instances of time in units of the energy 
density e. The direction of the gauge fields is in the transverse xy-plane and 
the momenta are parallel to the z— axis. This correlation function is typically 
associated to an occupation number distribution divided by frequency. 
Because of the oblate initial conditions the distribution along the 2;— axis 
is characterized by small values at early times. Exponential growth due to 
plasma instabilities leads to a rapid population of long-wavelength modes, 
while fluctuation effects induced by the growth in the lower momentum modes 
lead to amplification in a broad range of momenta. The exponential growth 
stops around t ~ (200 — 250)e~^/^, and we concentrate on the subsequent 
evolution. The behavior at earlier times is analyzed in detail in Ref. |4j. 

After the exponential growth period the dynamics slows down 
considerably. As time proceeds, more and more ultraviolet modes approach 
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Figure 2: Simulation results for transverse (square) and longitudinal 
momenta (circle) at a later time. The fit indicates a power-law behavior. 

a power-law described by 



with the "occupation number exponent" k. Fig. [2] shows the simulation 
results at time t = 4000e^^/'* for transverse (square) and longitudinal 
momenta (circle). Given that longitudinal and transverse momentum modes 
were populated very differently at initial time, the differences reflect the 
comparably small remaining anisotropy at this time. The solid line represents 
a power-law fit (jl]) with 



which matches the data in a significant momentum range rather well. 

In order to obtain the fit value for n from the simulations, we first 
employ a least-square fit to the transverse and longitudinal results separately 
at a given time. Subsequently, we average the results for the transverse 
and longitudinal momenta. For the longitudinal momenta we choose the 
fit range Pz € [0.6, 2] e^^^. For the transverse-momentum results, which 
display power-law behavior in a larger domain at earlier times, we use 
the range pj_ E [0.6, 3] Conver gence in the infrared proceeds rather 

slowly. However, changing the lower bound for the fit range from 0.6e^/^ 
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K = 1.35 ±0.02 , 



(5) 
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Figure 3: Fits to a power-law behavior with exponent k taken at different 
times. The dashed hnes represent perturbative Kolmogorov exponents. 

down to 0.2e^/^ leads only to changes in n which are comparable to the 
statistical error. Fig. [3] displays the fit values for function of time. 

The comparably large error for t < lOOOe"^/'^ comes from the insufficient 
isotropization at these times. Subsequently, the values for k change very 
little as a function of time and the errors become small. The dashed lines 
correspond to the results for perturbative Kolmogorov exponents for energy 
(5/3) and particle (4/3) cascades, which is explained in Sec. [31 The late- 
time behavior for t > 4500e~^/^ is consistent with a slow approach to a 
classical thermal equilibrium distribution with exponent one (solid line). A 
remarkable agreement between simulation results and k = 4/3 is observed in 
the quasi-stationary period around lOOOe"^/"^ ^ ^ ^ 45006^^/"^. We obtain the 
best fit value dSD by taking the time-average in this quasi-stationary period. 

We emphasize that the errors quoted are purely statistical. The different 
dynamics for longitudinal and transverse modes leads to the fact that the 
longitudinal distribution corresponds to somewhat larger estimates for k than 
what is inferred for transverse modes at not too late times. If this difference 
is employed to estimate the systematic error for k, the error will be less than 
ten per cent. 

Our results seem to be incompatible with the value two for k reported 
in Refs. [HI |18j based on Vlasov equations. Unfortunately, the origin of 
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that value is not fully understood so far. In particular, a perturbative 
analysis in terms of Kolmogorov wave turbulence similar to the one in 
Sec. [3] seems not to be able to explain it. There is still no consensus 
about the underlying reasons for the different numerical estimates of k 
using Vlasov and classical-statistical approaches, respectively. It has been 
pointed out that in the Vlasov treatment the hard modes represent static 
sources, which do not isotropize. In contrast, total energy is conserved in the 
classical-statistical lattice gauge theory simulation and all modes isotropize at 
sufficiently late times. However, for weak coupling and assuming a sufficiently 
large separation of scales both approaches are expected to agree. One could 
conclude that simulations on very much larger lattices to achieve the clear 
separation of scales underlying the Vlasov approach might be required to 
compare the two. We verified that simulations on lattices as big as 256^, with 
much less statistics, show no indication for quantitatively relevant changes 
of the K value. We found that already simulations on 64^ lattices accurately 
reproduce ([5]), which indicates that the results are rather robust. 

Clearly, none of these approximations are sufficient to quantitatively 
address the asymptotic late-time behavior, which should finally be 
characterized by a Bose-Einstein distribution for the gluons. While for 
low momenta the employed classical-statistical simulations are expected 
to give an accurate description for sufficiently high energy densities or 
occupation numbers, the high-momentum behavior will be altered by 
quantum corrections. In the following, we address aspects of this question 
using resummed perturbation theory at high momenta. 

3 Resummed perturbative quantum theory 

For quantum fields one can define two independent two-point 

correlation functions out of equilibrium, which may be associated to the 
anti-commutator and the commutator 

F^ix^y) = \{{Al{x),A'M}) , ptAx,y)=i{[Al{x),Al{y)]), (6) 

respectively. Loosely speaking, the spectral function p determines which 
states are available, while the statistical propagator F contains the 
information about how often a state is occupied. The spectral function 
is related to the retarded propagator G'(ij) and the advanced one as 
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Figure 4: Gluon part of the one- loop contribution to the self-energy with 
(2PI) resummed propagator lines. 



P = C^{R) ~ G(^A)- A tremendous simplification of thermal equilibrium is that 
the spectral and statistical functions are related by the fluctuation-dissipation 
relation, which is not assumed here [T6] . 

In the following we will show how a power-law behavior (jlj) can be 
explained for sufficiently high momenta in resummed perturbation theory, 
if quantum corrections are neglected. Since the scaling solution is time and 
space translation invariant, the correlators in ([6]) can be Fourier transformed 
to F{p) and p{p) with four-momentum p = (p°,p)EI Similar to ([6]) we 
consider statistical, 11(^)5 and spectral, n(p), components of self-energies 
defined as (TB] 

K)aM = G^i^lip) - G^ifM , (7) 

where summation over repeated Lorentz and color indices is implied. The 
translation invariant propagators ([6]) and self-energies obey the nontrivial 
identity [13] 

which can be directly verified by plugging in the above definitions. This 
equation is well-known in nonequilibrium physics and will be the starting 
point for our calculation. In the language of Boltzmann dynamics it states 
that "gain terms" equal "loss terms" for which stationarity is achieved |16]. 
Thermal equilibrium trivially solves ([8]), which we do not consider in the 
following. Instead, we will look for possible non-thermal stationary solutions 
in perturbation theory relevant at sufficiently high momenta. 

To exemplify the evaluation of the self-energies, we consider the one-loop 
diagram of Fig. H] appearing in the two-particle irreducible (2PI) effective 
action scheme, where self-energies are expressed in terms of self-consistently 



^We introduce a —i in Fourier transforms of the spectral (p-) and retarded/advanced 
components, such as p{p) = —i J d'^x e^P'^'^'^ p{x), while F{p) — f d^x e^''''^^'^ F{x) |T6] . 
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dressed propagators while vertices remain undressed. Using F^^ = S'^'^F^^, 
and p^^ = at one-loop order the gauge field contributions to the color- 

diagonal self-energies can be written as [19] 

nf;)(p) = -g' [ {27rr6^'\p + k + q)\Ff,MFs,{q) + lha{k)p5,{q) 

J kq L . 

nf;)(p) = g' [ {27r)U^'\p + k + q)\Ff,^{k)ps,{q)+Ppa{k)Fs,{q) 

J kq 

X V^^"^''^^^(p,fc,g), (9) 

with the notation = J d'^k/{2TT)'^. The Lorentz and momentum structure 
of the two three- vertices appearing in the one- loop expression for SU(2) gauge 
theory is contained in 

V^-^''^^\p,k,q) = [g^»{p-ky + g^^{k-qr + g^^{q-pr] 

X [g-^f^^p-kY + gf^^k-qy + g-'^q-pf]. (10) 

In accordance with (jl]), we are looking for scaling solutions which behave 

as 

F^,{sp) = \sy^^^^^F^,{p) , p^,{sp) = \s\-'sgn{s)p^,ip) (11) 

under rescaling with the real parameter s. The spectral function is 
antisymmetric, p^uip) = —Pufi{—p), while F^y{p) = Fy^{—p) and we assume 
symmetry in Lorentz indices. As in Ref. [7j, scaling solutions can be efficiently 
identified by integrating ([8]) over external spatial momentum p and suitable 
scaling transformations for coordinates. This is similar to what is done in 
the context of weak Kolmogorov wave turbulence for scalar models using a 
Boltzmann equation [5], and in this way the problem is reduced to simple 
algebraic conditions for exponents. Here, the analysis is complicated by the 
fact that for the perturbative result to order g"^ the two terms in ([8]) vanish 
separately on-shell, i.e. for p^ = ±|p|. The non- vanishing contributions of 
the diagram in Fig. |4] to order g^ can be conveniently analyzed by writing 
dZD as 

^..(p) = G'i^Hp)n(^")(p)G^^^(p) , P..(p)=G'(^nP)<^(p)Gji^(p)- (12) 

This is used to replace one propagator line in the one-loop graph of Fig. HI 
To order g^, the self-energies appearing in ( fT2l) can then be taken as (Q. 
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Following a well-trodden path [H [121 El E], we obtain by integrating over 
spatial momentum p from the stationarity condition 



= 9' 5^'\p + k + q)5^'\k + l + r)V^s^,,rxip,k,q)Vo.S'-y',(3r'y{k,l,r) 

J pqklr 



p. 



+ 



X 



pu 



p^ 



fP' 1 1 



p^ 

sgn I -7 I + 



pv 



pu 

sgn| 



sgn I — I - 1 



• (13) 



Here 

A = 4 - 3k (14) 

with A = A + 2k. In f ll3p the first term in brackets, ~ FF, vanishes for 
A = — 1, and for A = in the on-shell limit for F and p such that only 2 2 
scattering contributes. With f|T^ one observes that these correspond to the 
well-known Kolmogorov exponents for energy and for particle cascades, 

5 4 
- , or K = - , 
3 ' 3 ' 



K 



(15) 



respectively. In (1131) the second term in brackets, ~ pp, does not vanish with 
f fTSj) . Doing the equivalent calculation starting from the classical-statistical 
gauge theory instead of the quantum theory leads to the same stationarity 
equation ([T5]) . however, without this term ~ pp [TH]. Correspondingly, 
neglecting quantum corrections one observes rather good agreement of the 
perturbative solution k = 4/3 with the value ([5]) indicated by the full 
numerical solution of the classical-statistical gauge theory. A similar analysis 
can be performed for the two-loop gluon diagrams. 

For sufficiently high momenta, where characteristic occupation numbers 
are of order one, the quantum corrections appearing in fll3p become relevant, 
which preempts a scaling solution f|T5l) in the far UV. On the other 
hand, parametrically the perturbative estimate certainly breaks down for 
occupancies as large as Oll/g"^), where order-one corrections to the self- 
energies occur at any order in a loop expansion. As a consequence, one 
expects a window of momenta for which the scaling solution (fTHj) can describe 
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quantitatively the dynamics. In view of this, it is remarkable that the 
simulation results shown in Fig. [2] exhibit for low momenta only small 
deviations from the perturbative behavior. Given that the numerical results 
are rather insensitive to going to 256^ lattices, simulations on much larger 
lattices might be required to settle the important question of whether a new 
nonperturbative infrared fixed point exists in nonequilibrium QCD, similar 
to what is observed in scalar theories [7]. 

4 Conclusions 

Despite important differences the discussion of the thermalization process in 
heavy-ion collisions and cosmology after inflation shows remarkable similar- 
ities. In both cases nonequilibrium instabilities can lead to a fast period of 
exponential growth of occupation numbers. This is followed by a slow period, 
where the quantitative agreement concerns characteristic exponents or scal- 
ing functions. Our numerical as well as analytical results indicate that the 
nonabelian gauge theory can lead to perturbative Kolmogorov scaling expo- 
nents known from scalar models. In QCD simulations for large volumes an 
ambitious analysis in terms of gauge invariant quantities such as stress-tensor 
correlation functions would be desirable for a quantitative understanding of 
the nonperturbative infrared dynamics. One can hope that this leads to 
an improved theoretical understanding of the plasma's temporal evolution, 
which is crucial for the successful outcome of the experimental heavy-ion 
program. 
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program on " Nonequilibrium Dynamics in Particle Physics and Cosmology" 
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